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We determine the joint distribution of the exit times from a two- 
dimensional strip of a bivariate diffusion process. We consider two 
different situations; crossing or absorbing boundaries. In both cases, 
this distribution depends on the transition density of the process 
constrained to evolve under the boundary. Solving a two-dimensional 
Kolmogorov forward equation, we explicitly derive these two quan- 
tities for a bivariate Wiener process with drift and non-diagonal co- 
variance matrix. Explicit expressions for other diffusion processes in 
presence of either absorbing or crossing boundaries are not available. 
We propose a numerical algorithm, which is shown to be convergent. 
A comparison between theoretical and numerical results for Wiener 
and an illustration of the numerical approximation for a bivariate 
Ornstein-Uhlenbeck process in presence of absorbing boundaries are 
carried out. 

1. Introduction and motivation. The first passage time (FPT) problem for one-dimen- 
sional stochastic processes has been widely investigated through simulation, analytical and nu- 
merical methods [3, 7, 14, 17]. Besides its mathematical interest, the derivation of the FPT 
distribution is relevant in different fields, e.g. neuroscience, reliability theory, finance, and epi- 
demiology. In neuroscience, FPTs describe the times when the neuron releases an electrical 
impulse, called spike. In reliability theory, FPTs model the epochs when a crash of an object 
happens. In finance, FPTs describe the time when a bond or a stock reaches a certain value 
and it is profitable to sell or buy. In epidemiology, FPTs describe the times when an epidemic 
reaches a threshold level, causing a major disease outbreak. Connections between neurons, com- 
mon shocks and direct interaction between objects, dependencies between stocks in the same 
portfolio or belonging to the same market and interactions between populations suggest the 
presence of dependencies between FPTs. Therefore, it is of interest to extend the FPT problem 
to more general scenarios. 

Our aim is to solve the two-dimensional FPT problem of a bivariate diffusion process in 
presence of crossing or absorbing boundaries, motivated by problems arising in the framework 
of neuronal network modeling. This paper is inspired by [13], where the FPT of a two-dimensional 
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Brownian motion without drift, in presence of absorbing boundaries, was derived. The bivariate 
Wiener was proposed as an oversimplified model of a neural network. A more realistic model is 
proposed in [22] , where a multivariate Ornstein-Uhlenbeck (OU) process is introduced to describe 
the dynamics of k neurons. This model is obtained as a diffusion approximation of a multivariate 
Stein process. The presence of common inputs between neurons determines dependence in their 
membrane potential evolutions and in their spike times. A bivariate version of this model is 
investigated through simulations in [21]. 

For one-dimensional processes, transition densities are known for Wiener, OU and square root 
processes (also called Cox-Ingersol-Ross or Feller process), and all those cases where an explicit 
solution of the Kolmogorov forward equation is available [9, 18, 19]. The transition density of a 
process in presence of absorbing boundaries and the FPT density are generally unknown, except 
for Wiener [20] and a special case of the OU process [11]. As an alternative approach, numerical 
methods can be applied [20]. 

For multivariate processes, the transition density is known for Gaussian diffusion processes 
[4]. Furthermore, in [5], the FPT of one component of a bivariate Gaussian diffusion process 
through a constant boundary is studied. However, neither the transition density in presence of 
absorbing boundaries, nor the joint FPT density are available. Our aim is to compute them in 
the two cases of crossing or absorption at the boundaries. In presence of absorbing boundaries, 
we assume that the first component reaching its threshold, which we call the faster component, is 
absorbed there, while the other component, which we call the slower component, independently 
pursues its evolution till the epoch when it attains its boundary. Thus the considered process is 
a bivariate diffusion until the crossing of the faster component and then it becomes a univariate 
diffusion. For crossing boundaries, we assume that the faster component pursues its evolution 
after having reached its threshold, together with the slower component, which evolves until it 
crosses its boundary. 

In Section 2 we introduce notation and mathematical background used throughout the paper. 
In Section 3, we calculate the joint distribution of the FPTs of a bivariate diffusion process 
in presence of crossing or absorbing boundaries. In both cases, the distribution depends on 
the unknown joint probability densities of the slower component constrained to be below its 
threshold, and of the FPT of the faster component. We show that these unknown conditional 
densities solve a system of Volterra-Fredholm integral equations. In Section 4 we present a 
numerical algorithm to solve the system. The obtained numerical solution is then used to evaluate 
the joint distribution of the FPTs. In Section 5 we study the order of convergence of the error of 
the proposed algorithm. In Section 6 we determine explicit expressions of the joint distribution of 
the FPTs of a bivariate Wiener process with constant drifts and non-diagonal covariance matrix. 
In particular, the transition density of the bivariate Wiener process in presence of absorbing 
boundary conditions is explicitly derived as solution of a bivariate Kolmogorov forward equation. 
The results in Section 6 extend and correct those in [6, 12, 13]. Finally, in Section 7 we apply the 
numerical algorithm to evaluate the joint FPT density of bivariate Wiener and OU processes in 
presence of absorbing boundaries, comparing numerical and theoretical results for Wiener. 



FPT FOR BIVARIATE DIFFUSION PROCESSES 



3 



Throughout the paper, bivariate diffusion processes are considered. To apply the algorithm and 
numerically evaluate the joint FPT density, we are restricted to processes with a known bivariate 
transition density, e.g. Gaussian diffusion. 

2. Mathematical background. Consider a two-dimensional time homogeneous diffusion 
process X = {(X±, X2)(t);t > to}, originated at time to hi X(io) = xo = (#01, £02), solution of 
the stochastic differential equation 

(2.1) dX(t) = ix(X(t))dt + £ (X(t)) dW(t). 

Let B = (Bi,B 2 ) G R 2 be a two-dimensional boundary, with Bi > xoi,i = 1,2 and denote Tj 
the random variable FPT of Xi through the boundary Bi, defined by 

Ti = m£{t>to:Xi(t)>Bi} * = 1, 2. 

Before the first exit time from the strip (—00, B\) x (—00,-82), the process X evolves under the 
boundary B, and we denote it as X a , i.e. 

X a = {X(t);te [t ,mm(Ti,T 2 )}} . 

Similarly, we denote Xf the unidimensional process Xi evolving under the boundary Bi, before 
time Ti, i.e. 

X? = {Xi(t);t€ [t ,Ti\}, i = l,2. 

If Z is a /c-dimensional process, we denote its transition probability by Fz(x, t|y,s) = P(Z(£) < 
x|Z(s) = y), its survival function by Fz(x, t\y, s) and the transition probability density function 
(pdf) by /z( x , *|y , s) for s < t,x, y G Throughout the paper, Z = X or Z = X a , for k = 1 
or 2. 

For s < t, we denote fx a \x a { x ii A x jit\ Yj s ) the conditional pdf of Xf given Xj, for i,j = 1, 2,i 7^ 
j, defined by 

fx-\x-{x i ,t\x j ,t;y,s)dx i = F(Xf(t) G dXi\Xf(t) = ^,X a ( S ) = y), 

and fx?\Tj { x i\t'j Vii s ) the transition pdf of Xf conditioned on Tj when Xf starts in j/j at time 
s < Tj, defined by 

(2.2) f x? \ T . ( Xi \t; y t , s) d Xi = P (Xf (Tj) G dxi\Tj = t, Xf (s) = Vi ) . 

The joint pdf of Tj and Xf(Tj) when X a starts in y at time s < Tj is denoted by f(x°-,Tj)( x i-> t\Vi s ) 
and defined by 

(2.3) / W>T .) (si, a) dxidi = P (Xf (Tj) G dx t ,T 3 G <ii|X a (a) = y) . 
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We denote /(Xj,T) (( x i> *l(2/i> 2/j')> s ) the joint pdf of Tj and Xj(Tj), when X,- starts in yj, Xi 
starts in yi at time s, which is defined by 

f(x t ,Tj) ((xi, Bj),t\y, s) dxidt = F (Xi (Tj) £ dxi, Tj £ dt\X (s) = y) . 

The pdf of T is defined by gx { (t\yi, s) = P(T £ dt\Xi(s) = yi). Our aim is to determine the joint 
distribution function of T = (T\, T2) for a process X starting in y < B at time s in presence of 
either crossing or absorbing boundary B. 

To simplify the notation, throughout we omit the starting position when y = Xo, and the 
starting time s when s = to- 

3. Joint distribution of (T±, T 2 ). In presence of absorbing boundaries, the joint distribu- 
tion of (T\,T 2 ) can be expressed in terms of the marginal FPTs densities and of the joint pdfs 
(2.3) as follows 

Theorem 3.1. Let X be a two-dimensional diffusion process with X(io) = x o an d B be a 
two-dimensional absorbing boundary with B\ > xqi and B2 > xq2- The joint distribution of 
(Xi,T 2 ) is 



(3.1) 



P(Ti <t u T 2 <t 2 ) 



^2 [ [ ([ 9T 3 (sj\xj,Si)ds^\ f (X a Ti) (xj,Si)dxjdsi 

-I.A^-^Jto J -CO \Jsi / 



Proof. 



P(Ti <h,T 2 <t 2 ) 



2 . t . 

= / P(^l < ti, T 2 < t 2 \Ti < Tj, Ti = Si )¥(Ti e d Si ,Ti < Tj). 

Conditioning on the value of the component which has not yet reached its boundary, at the time 
when the other component crosses its boundary, we get 

P(Ti <h,T 2 <t 2 ) 
2 r u nBj 



Y F(Tj<tj\T i =s i ,X^(s i )=x j )F(X^(s i ) £ dxjlT^sjFlTi £ dsi,Ti<T. 

2 r-ti pBj 

Y / / HTj < tj\Xf(8i) = Xj )P(Xf( Si ) £ dxj,Ti £ dsi), 



where the last equality holds because X and thus X a are Markov processes. □ 
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Remark 3.1. The expression 

(3.2) f( X ?,Ti)(xj>Si)dxjdsi = fx^\Ti{xj\si)g Ti {si)dxjdsi 
can be plugged into (3.1), obtaining 

P(Ti < h,T 2 < t 2 ) 

(3.3) = S / / (/ 9T j (aj\xj,8 i )daA fxa^Xjls^grMdxjdai. 

This expression is useful when gx t is known, because it allows to rewrite (3.1) in terms of the 
unknown function fx a \T- 

From Theorem 3.1, it follows 

Corollary 3.1. The joint density of (Ti,T 2 ) for a two-dimensional diffusion process in 
presence of an absorbing boundary B is given by 

(3.4) P(Ti G dh,T 2 e dt 2 ) = / gr^tjlxj^^f^x^T^ixj^^dxjdtidtj 

-:i L/ •• ,j x 

In presence of crossing boundaries, the joint FPT distribution of a diffusion process is given 
by the following 

Theorem 3.2. Let X be a two-dimensional diffusion process with X(io) = xo and B be a 
two-dimensional crossing boundary with B\ > xoi and B 2 > xq 2 . The joint distribution ofT is 

/ P(Ti <h,T 2 <t 2 ) 

= / / fix^T^XuB^^s^Bi.x^.s^f^xa^x^s^dxidsjdxjdsi. 

ij — l-i-t-j-lto J -CO J Si J -co 

Proof. Before the first crossing time from the strip, the bivariate diffusion process behaves 
as the bivariate diffusion process. Therefore, from Theorem (3.1), it holds 
(3.6) 

P(Ti < h, T 2 <t 2 )= / / P ( T i < h\ T i = *i (*<) = Xj )F(Xf(si) E d Xjl Ti E d Si ). 

• ■ i.j / j Jtn J —CO 



After the first exit time from the strip, the faster component i pursues its evolution, starting in 
Bi at time Si, while the component j starts in Xj(si) and has its FPT at time Tj. Therefore, we 



6 L. SACERDOTE, M. TAMBORRINO AND C. ZUCCA 

have 

F{Tj < tj\Ti = Si,X?( Si ) = Xj ) = I ' F(Tj G d Sj \Ti = Si ,X^{ Si ) = Xj ) 

J Si 

f'tj /*oo 

= / / F(Tj € d8j,Xi(8j) £ Xi\Ti = Si,X°(si) = Xj ) 

J Si J — oo 
f'tj /*oo 

(3.7) = / / F(Tj £ d Sj ,Xi( Sj ) £ Xi\Xi{si) = B u X«(si) = Xj ) 

J Si J—oo 

where the last equality holds because X is a Markov process. The theorem follows from plugging 
(3.7) into (3.6). □ 

The density f(Xi,Tj) m (3-5) is not explicitly known, but it can be numerically evaluated as 
done in [5], provided that the bivariate transition density is known, e.g. for bivariate Gaussian 
diffusions. Note that the joint FPT distributions (3.1) and (3.5) of a bivariate diffusion process 
in presence of absorbing and crossing boundaries, respectively, depend on f(x a ,T t )- I n Section 
6, we explicitly derive f(x a ,Ti) f° r a bivariate Wiener. For other processes, the unknown density 
f(x a ,T x ) solves a system of Volterra-Fredholm integral equations [1]: 

Theorem 3.3. Let X be a bivariate diffusion process with X(to) = xo and let B be a two- 
dimensional boundary with B\ > xoi and B2 > xo2- The joint transition pdfs f(x a ,Tj)> f or 
i,j = 1,2; i 7^ j, are solutions of the following system of Volterra-Fredholm first kind integral 
equations 

ft r B 2 



Fx((x 1 ,B 2 ),t)= j j 2 Fx{(x 1 ,B 2 ),t\(B 1 ,y),T)f {x * tTl) (y,r) dydr 

J tg J — oo 
ft fB 1 

(3.8a) +/ / Fx({ Xl ,B 2 ),t\(y,B 2 ),T)]f {x « >T2) (y,T)dydT; 

J t() J — oo 

F x ((B 1 ,x 2 ),t)= f f 2 Fx({B 1 ,x 2 ),t\(B 1 ,y),T)f {X a,T 1 )(y,T)dydT 

J to J—oo 
ft fB x 

(3.8b) +/ / Fx((B u x 2 ),t\(y,B 2 ),T)f {xtjT2) (y,T)dydT, 

J to J — oo 

where x\ > B\ and x 2 > B 2 . 

Proof. Let us consider the exit times of the process X. The survival distribution of X, for 
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xi > B\ and X2 > B2, is given by 

Fx(x,t) = P(X(i) >x,Ti <T 2 )+P(X(t) >x,Tx > T 2 ) 

= / [ 2 F(X(t)>x,T 1 <T 2 \T 1 =T,X 2 (T 1 ) = y) 
Jto J — OO 

• P(X 2 (ri)€dy,Ti€dr,ri<r 2 ) 
+ / /" 1 p(X(i)>x,T 1 >T 2 |r 2 =T,X 1 (T 2 ) = y) 
■ F(X 1 (T 2 )edy,T 2 edT,T 1 >T 2 ) 

= ! ! 2 ¥(X(t)> X \X(T) = (B 1 ,y))f iX a tTi) (y,r)dydT 
Jto J — OO 

(3.9) + /" P (X(t) > x|X(r) = (y, B 2 )) / ( x f ,T 2 ) fo, r) dycfr, 

where the last equality follows from the strong Markov property. Then, the theorem (3.8) follows 
by choosing x = (xi, B 2 ) and x = (B\,x 2 ), respectively. □ 

Corollary 3.2. Differentiating (3.9) with respect to x and since 
/to f-L /(X?,T 2 )(2/> r)dydT = 1, it follows 

/x(x,t) = / / fx(^,t\(B 1 ,y),r) f {X a Ti) (y, T )dydT 
J to J —00 
rt rBi 

+ 11 /x(x,t|(y,S 2 ),r)/(x»,r 2 ) {y,T)dydr. 
Jto J — 00 

Remark 3.2. Assume to have a bivariate process with a component which is reset to its 
starting value whenever it crosses its boundary, and then it restarts. Meanwhile, the other com- 
ponent pursues its evolution. The intertimes between two consecutive FPTs of the same com- 
ponent are independent and identically distributed and thus each component is marginally mod- 
eled by a renewal process. The joint FPT distribution depends on f(x<?,Ti) an d f(Xi,Tj)- The 
second density describes the evolution of the process in (Ti,Tj) when T{ < Tj: X starts in 
X(Ti) = (^Xj(Ti), XQ^j , where the crossing component Xi is reset to its initial value XQi. Then 
the slower component j has its FPT at time Tj, while Xi evolves, possibly with further crossings 
of Bi (with corresponding reset) before Tj, with i,j = 1,2. 

4. Numerical method and its convergence property. For a diffusion process, the den- 
sity /(Xf,Ti) is generally unknown. Therefore the system (3.8) cannot be analytically solved 
and neither the joint FPT distribution (3.1) nor (3.5) can be explicitly derived. To solve (3.8) 
and therefore approximate /(x a ,T,)> we propose a numerical method which does not depend on 
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whether the boundaries are crossing or absorbing. Consider an assigned two-dimensional time 
interval [0, 0i] x [0,@2], with ©i,©2 € M + . For each component i = 1,2, let hi and be the 
time and space discretization steps, respectively. On {[— oo,-E?i] x [—00,^2] x [0, @i] x [0,02]} 
we introduce the partition {(y Ul , y U2 )',t kl , tk 2 } where t ki = hhi is the time discretization and 
y Ui = Bi — UiTi is the space discretization for fcj = 0, . . . , iVj, Nihi = 0j, Ui G N, and i = 1,2. 
To simplify the notation, we consider hi = h 2 = h and ©i = ©2 = 0, implying N\ = N 2 = N, 
k\ = k<i = k and thus t kl = t k2 = t k , for k = 0, . . . , N. 

Let f(x^,T 2 )(y^j) denote the approximation of f(x a ,T 2 )(y>tj) due to the time discretization 
procedure. We approximate the time integrals of (3.8) through the Euler method [15], obtaining 
a system of integral equations for f(x a ,T 3 ){Vi t), i = 1,2. For x\ < B\ and X2 < B2, we get 



Fx((x 1 ,B 2 ),t k )=h^2 / F ^(i.xi,B 2 ),t k \(Bi,y),tp)f (X a >Tl) (y,tp)dy 





k r B 1 



(4.1a) 



+ h '^Z I F ^(( x i^ B 2),tk\(y,B2),t p )f {X a T2) (y,t p )dy; 



Fx((B 1 ,x 2 ),t k )=h^2 / Fx((B 1 ,x 2 ) ,t k \(B 1 ,y),t p )f( X a Ti) (y,t p )dy 





k r B 1 



(4.1b) 



+ h ^2 / F ^(( B i^ x 2), t k\(y,B2),tp)f {X a T2) (y J t p )dy. 



Let 1a denote the indicator function of the set A. Then 



(4.2) 



Fx{(B 1 ,x 2 ),t k \{B 1 ,y),t k ) = t{y >X2 y, 
F K ((B 1 ,x 2 ),t k \(y,B 2 ),t k ) = 0; 
Fx({x 1 ,B 2 ),t k \{B 1 ,y),t k ) = 0; 
Fjc((xx,B 2 ),t k \(y,B 2 ),t k ) =; l{ y>xi y 
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Plugging (4.2) into (4.1) and differentiating with respect to Xj,j = 1,2, we get the system 
dF^((x u B 2 ),t k ) 



Q hf {X a T2) ( Xl ,t k ) 

" Ba dF^((x 1 ,B 2 ),t k \(B 1 ,y),t p 



— r^ dF^x 1 ,B 2 ),t k (B 1 ,y),t p ) f 

+ h l^J_ oo ox~ x fasau (v. *p) d y 



p=0 

— 1 ^ aFx((xi,i?2),^|(y,i?2),^ 



( 4 -3a) + h l^ «i /(x?,t 2 ) (y, ip) dy; 

p=Q J-oo 1 

0F x ((i?i^ 2 ),t fc ) . 

^ [ B * dFx((B 1 ,x 2 ),t k \(B 1 ,y),t p ) f 

+ h 2^ ^ /(xlto (y. t P ) dy 

(4.3b) + h 2^ ^ f(x?,T z ){y,tp)dy. 

p=0 J-co <Jx 2 

Discretizing the spatial integral and truncating the corresponding series with a finite sum, we 
obtain 

dF x (( Xl ,B 2 ),t k ) ~ 

0^- = nj(x?,T 2 ) {xi,t k ) 

^ ^ dFx((xi,B 2 ),tk\(B 1 ,y U!i ),tp) ~ 
+ hr ^2^ qT hx$,Ti) (yu 2 ,t p ) 

p=0 U2=0 

( a a \ , h dFxdx!, B 2 ),t k \(y Ul , B 2 ),t p ) - 

p=0ui=0 1 
^ «/(X?,2i) (X 2 ,*fc) 

fc__i ^ aFx((-Bi,x 2 ),tfc|(-Bi,j/ M2 ),t p ) - 

p=0« 2 =0 2 

(4.4b) + /in 2^ 2^ fli f(xi,T*){yu X ,tp)- 

p=0« 1= WX2 

Here f(x a ,Tj) (j/>*) denotes the approximation of /(x<%t) (y>*) due to the time and space dis- 
cretization procedures and to the truncation of the infinite sums of the space discretization. 

Since f^ X a ,T j )(yu i , to) = 0, we set f(Xi,Tj)(yui,to) = 0. The following algorithm can be used to 
approximate f{Xf,Tj) in tne knots {(Vm > ^2); 
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1 _d_ 
h dx± 

1 d 



/(Xf.Ti) (yu 2 ,h) = - —Fx((B 1 ,x 2 ),t 1 ] 



xi=y m 



X2=yu 2 



Step k > 2 



/(Xf,T 2 ) (Uuntk) 



1 _3_ 



F x ((xi ) S 2 ) ) t*) 



xi=y u 



fc— 1 m.2 



p=0 i; 2 =0 
fc— 1 mi 



r iJ3 S /(A?,T 2 )(ifai>*p) [F x ((x!, B 2 ),t k \(y Vl , B 2 ),t p ))] 

p=0 Vl =0 1 



Xl=Hu 



1 _3_ 



Fx((B 1 ,x 2 ),t k ) 



X2=Vu 2 



fc— 1 m 2 



r2^2 ^2 h^,Ti) (yv 2 ,tp) Q^-Fx((Bi,x 2 ),t k \(B 1: y V2 ),t 



p=0 1>2=0 

fc— 1 mi 



3 



22=£/u 2 



- n. /(xf,T 2 ) ^F x ((-Bi,x 2 ),ifc|(?/ wl ,S2),ip)) 

Note that at time ti, /pq>,T 3 -) (y 5 *l) = f(x?,Tj) {VinM) in each knot 2/«r 

Remark 4.1. VFe choose the Euler method because it simplifies the notation and is easy to 
implement. More efficient schemes, e.g. trapezoidal formula, can be similarly applied, improving 
the rate of convergence error of the proposed algorithm. 

5. Convergence of the algorithm. Let E <yt \y Ui ,tk) denote the error of the proposed 
algorithm evaluated in the mesh points (y Ui ,tk), for k = 0, . . . , N, itj = 0, 1, . . . , rm,i = 1, 2. It 
is defined by 

(5.1) E^(y Ui ,t k ) = /(x?,Tj) (y«i>4) - /(X?,T.,-) (2/^,4) , hJ = l,2,i ^ j. 

Mimicking the analysis of the error in [8], we rewrite the error (5.1) as a sum of two errors. The 
first is given by e^\y Ui ) = f(x?,Tj) {y Ui ,tk)- f{Xf,Tj) {Vu z ,tk) and is due to the time discretization. 
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The second is given by E^ u . = f(x?,Tj) bju t ,tk)- f\x^T 3 ) {y Ui ,tk) and is determined by the spatial 
discretization and by the truncation introduced at steps k > 2. We start computing E^ u , through 
the following 



Lemma 5.1. It holds 



E 



(i) 



k-l 



Bi 



Ki, k ,p{{y ui , B 2 ), (y, B 2 ))f( X <?,T 2 )(y, t p )dy 



p=0 

- / K ltk;P ((y ul ,B 2 ),(B 1 ,y))f( X a :Tl) (y,t p )dy 

J —c 



-oo 
mi 



(5.2a) 



4 3 =y 

k,u 2 / j 
p=0 

rB 2 



v 2 =0 
k-l r r B, 



(5.2b) 



where the kernels are 



+ n^2 K 1Ap ((y Ul ,B2),(y Vl ,B 2 ))f(x?,T 2 )(yv 1 ,tp) 

m 2 

+r 2 K^Kpdyu^B^^B^y^))];^,^) (y V2 ,t p ) 

K 2 ,k,p((B 1 ,y U2 ),(y,B 2 ))f {X a >T2) (y,t p )dy 

K 2 ^p{( B i,yu. 2 ),{B 1 ,y))f {X a tTi) (y,t p )dy 
i 

+ ri K 2,k,p{( B ^V^)^ (jfoi> B 2 ))f( X a^ 2) {y Vl ,t p ) 

di=0 
m 2 

+r 2 ^2 K 2)kiP ((Bi,y U2 ),(Bi,y V2 ))f iX a Tl) (y V2 ,t f 



— CO 

mi 



v 2 =0 







(5.3) 



Ki, k , t ((y ui , b), (c, d)) = [F x ((xi, 6), * fc |(c, d), t) - F x ((*i, 6),t fc _i|(c, d), t)] 



K 2 , k , t {(a,y U2 ), (c,d)) = — [F^((a,x 2 ),t k \(c,d),t) - F^((a,x 2 ),t k - 1 \(c,d),t) 



xi=y ui 



x 2 =yu 2 



When t = t p , we write Ki jkjP instead of Ki^ kttp to simplify the notation. Here a, c G (— oo,i?i) 
and 6, de (—00,-62)- 
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Proof. Subtracting (4.4) from (4.3), we obtain 



k-i r- . 

f(Xf,T 2 ) ( x l^k) - f(X?,T 2 ) ( x l>tk) = ^ - / 

p=0 L J ~ 



dFx((x 1 ,B 2 ),t k \(y, B 2 ),t p ) f 

f(x?,T 2 ) {y,tp)dy 



B > dFx((x 1 ,B 2 ),t k \(B 1 ,y),t p ) 



dxi 



dxi 



/(X 2 «,Ti) (y,t p )dy 



Ml=0 



"" dFx((x 1 ,B 2 ),t k \(y Ul ,B 2 ),t p ) - 

J(Xf,T 2 ) {Vu^tp) 



dx\ 



(5.4a) 



+ r 2 2^ /(Xf I r 1 )(l/«a,*p) 

«2=0 



fe-1 



/(Xf,Ti) (x 2 ,tk) - /(Xf,Ti) (x 2 ,tk) = 52 ~ I 



Bl 5Fx(( J Bi,x 2 ) ) tifcK2/, J B2),tp) 



Sa dFx((B 1} x 2 ) ,t k \{ B uy),tp 



— oo 

A; mi 



<9x 2 

-/(X^.Tj) {y,tp)dy 



f{x?,T 2 )(y,t p )dy 



•sr^ dF-x((Bi,x 2 ),tk\(y Ul , B 2 ),t p ) , 
+ r iL &^ /(x?,T 2 )(2Au,i P ) 



p=0m=0 

/«2 



(5.4b) 



^ aF x (( J Bi,^2),tfc|( J Bi,y U2 ),t p ) 
+ r 2 2^ &^ f{X^T 1 )(yu2,tp) 



u 2 =0 
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If we rewrite (5.4) for k — 1, without making the (k — l)th term explicit, we obtain 

nBl dFx(( Xl ,B 2 ),t k \(y,B 2 ),t^ - 
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fc-i 

E 

p=0 



B2 dFx((x l ,B 2 ),t k \(B 1 ,y),t p ) 



~f{Xf,T 2 ) {y,t p )dy 



/(X^,Ti) {y,t p )dy 



dFx((xi,B 2 ),t k \(y Ul ,B 2 ),tp) ~ 

J{X<t,T 2 ) {Vu! , t p ) 



M1=0 
1112 



dx\ 



(5.5a) 



^ dF x ((x 1 ,B 2 ),t k \(B 1 

+ r 2 2^ TtZ f{X-,T 1 ){yu2,tp) 



u 2 =0 
k-1 r ,B 



E 

p=0 



dx\ 

1 dFx((B l ,x 2 ),t k \(y,B 2 ),t p ) 



0; 



dx 2 



B2 dF x ((B l ,x 2 ),t k \(Bi,y),t p ) 



f(x?,T 2 )(y,t P )dy 



oo 
k mi 



dx 2 



f(xs,Ti) (y^p) d y 



sr^ dF x ((B 1 ,x 2 ),t k \(y Ul ,B 2 ),tp) ~ , 

+ r iL ^ f{x<t,T 2 ){y^t p ) 



p=0 Ul =0 

1112 



(5.5b) 



5F x (( J Bi,x 2 ),tfc|( J Bi,y U2 ),t p/ 
+ r 2 2^ &^ f{X-,T 1 ){yu2,tp) 



u 2 =0 



0, 



due to conditions (4.2). Then, the result follows by subtracting (5.5) from (5.4) and setting 
Xi = y Ui , for i = 1,2. □ 

To prove the convergence of the proposed algorithm, we need the following conditions: 

Condition 5.1. For i, j = 1, 2, i ^ j, k = 1, . . . , N, p = 0, . . . , k — l,in = 0, . . . , m*: 

(i) K ijktP ((a,b),(c,y))f(xf,Ti) (y>tp) and K iAp ((a,b), (y, d))/(x?,T i )(y>*p) are ultimately mono- 
tonic in y; 

(ii) f(x a ,Tj){y-,tp) is bounded, belongs to L 1 and there exist positive functions C^xiy) € L 1 , 
Ci <2 (y) G L 1 , with x G (— oo,Bi) and y G (—00,-82), swc/t i/mt 

|-£Q,fc,p((a, &), (2/, d))\ < hC it i(y); 
\Ki t k, p ((a, b), (c, y))\ < hC ij2 (y), 



and Cj,i(0) and Ci >2 (0) are bounded; 
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( Hi ) for 1 = 1,2 
(5.6) 



I- 



Bi—mi{n)ri 



QAy) fix?,?,) {y^p) 



dy < iiLiTi 



as Ti — > and rrii{ri)ri — > oo, where iftn are positive constants; 
(iv) for I = 1,2, there exist constants Qi t i such that 



(5.7) 



Si 



0_ 
dt 



dy 



< Ql,i', 



(v) for I = 1,2, zi = (y Ul ,B 2 ) and z 2 = (Bi,y U2 ), 




Fx(zi,tp\(B 1 ,y),t)f\ x %,Ti)(y, t ) \t=r £ L 1 in y £ (-oo, B 2 ); 
Fx(zi,t p \(y,B 2 ),t)f( X a )T 2)(y,t) \ t=T G L 1 in y G (-oo, B x ); 

Fx{iht p \(Bi,y),t)f( X a Tl) (y,t) \ t=T G L 1 in y G (-oo,-B 2 ); 

F^{zi,t p \(y,B 2 ),t)f^ i T 2 )(y, t ) \t=r G L 1 in y e (-oo, Bt). 

The following theorem gives the convergence of the proposed algorithm. 

Theorem 5.1. If conditions 5.1 are satisfied, then 
(5.8) £«(y Ui ,^) = 0(/i) + 0(r), 

where r = max(ri, r 2 ). 



Of 
0_ 

dt 
d d 

dy U[ dt 

d d 

dy Ul dt 



(i) 

Proof. At first, we study the error E) ' due to the spatial discretization. It can be decom- 
posed as 



(5-9) 4i= A ^- B Z^ k - ] V 

Here, has the same expression of E^ u , in (5.2), replacing f(x a ,Tj) (y> tj) with f(x a ,Tj) (Vi tj)- 
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(i) 

Moreover, B y k ' u . is defined by 



(5.10a) 



(5.10b) 



k-l 

b£ ] =y 

K,Ul / j 

p=0 
rti2 



mi 



i) 



di=0 



+r 2 E K 1Ap ((y Ul ,B 2 ),(B 1 ,y V2 ))Eg 



2) 

V2 



v 2 =0 



B 



(2) 
k,U2 



k-l 



mi 



= E ri E ^W^i,^), (y^.Ba))^ 

p=0 L «i=0 

+r 2 E K 2 , fclP ((Bi,y U2 ),(Bi,y, 2 ))4% 

1)2=0 



The term A^, accounts for the approximation of the spatial integrals with finite sums. Hence 
we can split it into a first term A^'^ accounting for the discretization procedure and a second 
A^if) for the truncation of the series. 
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By definition of , we have 



14 



Ul I 



Ki t k, P ((yui , B 2 ), {y, B 2 ))/(x« ,t 2 ) (y, t P )dy 



(5.11a) 



U (2,a), 



E 

oo 

-ri ^ K ltkiP ((y Ul , B 2 ), (y Vl , B 2 ))f(x<?,T 2 )(yv 1 ,tp) 

t)i=0 

+ / Ki :kiP ((y Ul ,B 2 ),{B 1 ,y))f( X a Ti) (y,t p )dy 

_J — oo 

oo 

-r 2 ^ K hkiP ((y Ul ,B 2 ), (B 1 , y V2 ))f(x^,T 1 ) (yv 2 ,t p 

v 2 =0 
k-l 



E 

p=0 



K 2<k , P ((Bi,y U2 ), (y> B 2))f{x°,T 2 )(y, t P )dy 



+ 



r\ ^ K 2,k,p(( B l,yu 2 ),(yv 1 ,B 2 ))f^ X ^ T2) (y Vl ,t p ) 
di=0 

^2,fc, P ((-Bi,yM 2 ), (B 1 ,y))f (X a >Ti) {y,t p )dy 



-oo 
oo 



(5.11b) 



-r 2 ^ K 2 ,k, P {{.B 1 ,y U2 ), (B 1 ,y V2 ))f ( ^ X a Tl) (yv 2 ,t p 

D2=0 



Let us focus on the terms in the first square brackets in (5.11a). It holds 

Ki^ p {{y Ul , £2), (y, B 2 ))f {X a yT2) (y, t p )dy 



r i ^2 K i,kA(yui>B 2 ), {y Vl , B 2 ))f(x?,T 2 )(yv 1 ,t l 

di=0 

Ki, k ,p{{y Ul , B 2 ), (y, B 2 ))f( X ?,T 2 )(y, tp)dy 



< 



Bi 



(5.12) 



< /l 

'Bi-n 



Ci,i(y) f(xf,T 2 )(y,t P ) 



dy 



where we used condition 5.1(i) and eq. (3.4.5) in [10] in the first inequality and condition 5.1(h) in 
the second. Note that the numerical approximations f(xf,Ti)(y^ t p ) can be rewritten as a function 
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of f(x a ,T i ){Vi'tp) an d Ki^, p {a,b,c,d), as shown in Remark 2 in [8]. Then, thanks to condition 

5.1(h), it follows that f{x a ,Ti){y^p) ls bounded. Moreover, the integrable function C\^{y) on 

the compact interval [B\ —r\,B\\ is bounded. Thus Ci ) i(y)\\f{x a ,T 2 )(y^p)\ — Vi,l f° r a positive 
constant rji, which yields (5.12). A similar procedure can be done for the terms in the second 
square brackets in (5.11a) and for those in (5.11b), obtaining 



(5.13a) 
(5.13b) 



fc-i 

I^Tj I < (nvi,i + r 2 r?i, 2 ) ^ h = (rir/1,1 + r 2 r/i )2 )tfe-i; 

fc-i 

l^l 2 ^ I < + ^2^2,2) ^2 h = ( ri? ? 2 - 1 + r 2V2,2)tk-i- 

p=0 



Here 77^ are positive constants given by 

Cl,i(y)f(X?,T 3 )(y,tp) < 

for i,j, I = l,2,i ^ j. Let us now consider the error At^ . Using condition 5.1 (i) , eq. (3.4.5) in 
[10] and then conditions 5.1(h), 5.1 (hi) in sequence, we get 



\A 



(1,6)1 

k,ui I 



fc-1 

E 

p=0 



n K hkA(yui, B 2),(yv 1 ,B 2 ))f(x?,T 2 )(yv 1 ,tp 



i)i=mi+l 



+r 2 ^2 K l,kA(yui^ B 2),(Bi,y V2 ))f iX a T ^(y V2 ,t p 

V2=m2 + 1 

k-l 



< 



+ 



p=0 



Bi—miri 



Ci,i{y)f(x?,T 2 )(y,tp)dy 



B2—m2T2 



Ci, 2 (y)/(xf,Ti) (y^p) d y 



(5.14a) 
(5.14b) 



< (^1,1^1 + ^1,2^2)^; 
\A%®\ < (^ 2l iri+V2, 2 r 2 )t fc , 



where (5.14b) is obtained as (5.14a). 

Prom (5.13), (5.14) and r = max(ri, r 2 ), we get l^-j^J < rGitk, where Gi, i = 1, 2 are positive 
suitable constants. Using these bounds in (5.9) and observing that B^' u , in (5.10) involves the 
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(i) 

errors E y p! Vi for < p < k - 1, we get a system of inequalities 



\E { ^ 1 \<G 1 rt k + rJ2 



k-l 



p=0 



mi 



Y,\Ki,kA(yu 1 ,B 2 ),(y Vl ,B 2 ))\\E^ 1 \ 



vi=0 



in-2 



(5.15a) 



+ E \K hk , p ((yu 1 ,B 2 ),(B 1 ,y V2 ))\\E { S 2 



V2=0 



\E^ 2 \<G 2 rt k + rJ2 



k-l 



p=0 



nil 



Y,\K2,kA B i>yu*),{yv^B 2 ))\\EM x \ 



.^1=0 



in 2 



(5.15b) 



+ ^\K 2 ^(B u y U2 ),(B u y V2 ))\\E { p 2 l 

V2=0 



t'2 



We extend the method proposed in [8] to the system (5.15), that we solve iteratively as follows: 



l<J = 0:=r^; 
\E[%\<G t r tl =:rpf- 



(0. 



\E^l\<G ir t 2 



r Pi ] ^ \ K ^k, P {(yui,B 2 ) : (y vl ,B 2 ))\ 



Till 



di=0 



+rpi 2) Y, \K l!k>p ((y ul ,B 2 ),(B 1 ,y V2 ))\ 



v 2 =0 



(5.16) 



< r \dt 2 + TftpP + rPbi 



(2) 



\E%U<t G 2 t 2 + rffp\ x > + r%p\ 



2JX) , ^ R 2 (2) 



(1) 

=: rp\ 

(2) 

:= rp x 2 , 



where (5.16) holds due to condition 5.1(h) and eq. (3.4.5) in [10]. Here f3j,i,l = 1,2 are suitable 
constants, which do not depend on r and h. Iterating this procedure, (5.15) becomes 



(5.17a) 



Kl\<r 



k-l 



Git k + rY J {PlP p 1) +^ 

p=0 



(5.17b) \E ( k %\<r 



k-i 

G 2 t k + 

p=0 



(i) 

:= rp\ '\ 



(2) 

:= rp\ . 



Since t k < G, from (5.17) it follows 



k-l 



pf <G i @^rY,{^P p 1) +P\P l 

p=0 



P I ' 



i = l,2. 
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Then, by eq. (7.18) in [15], we get pf < GiQexp[{p\ + $)rt k ). Therefore 

(5.18a) < rGiSexp [(/3j + $)rt k ] ; 

(5.18b) l4 2 ij < rG 2 eexp + /3 2 2 )ri fe ] , 



19 



ht k r-Bi d 



implying \E^ U .\ = O(r). 

Now, we focus on the time discretization error e k ^(y Ui ). The error formulas for the Euler 
method are 

(5.19) *i,i, fc (fc) 

<*2,l,fc(fr) 



I-L m F x((ym,B 2 ) 1 t k \(y,B 2 ),t)f {X a :T2) (y,t)dy 



ht k rB 2 d 



JZ> WtMiVu, , B 2 ),t k \(B 1 ,y),t)f {X a Tl) (y, t) dy 



ht k rBi d 



/_« ( (Si , y U2 ), t fc | (y, fl 2 ), 0/(x f ,r 2) (y , *) dy 



t=T 



t=T 



t=T 



ht k rB 2 d 



mM(Bi,y U2 ),t k \(Bi,y),t)f {X a Tl) (y, t) dy 



where r 6 (0,0) and we used = /i£;. Rewriting (4.1) with the corresponding residuals and 
evaluating it in Xi = y Ui , i = 1, 2, we get 



(5.20a) 



FxiiVm ,B 2 ),t k \(y,B 2 ), t p )f( X <?,T 2 ) (y,t p ) dy 

P =v -°° 

+ h^2 Fx({yu 1 ,B 2 ),t k \(B 1 ,y),t p )f(x$,T 1 ){y,tp)dy 
p=o j -°° 

+ h,i,k(h) + 5i )2 ,fe(/i); 

k r-B 1 

Fyi((B 1 ,y U2 ),t k ) = / F x(( B i,y U2 ),t k \(y, B 2 ),t p )f( X a tT2 )(y,t p )dy 

P =o J -°° 
k ,.b 2 



+ h ^2 F x((Bi,y U2 ) ,t k \(B 1 ,y),t p )f( X a Ti) (y,t p )dy 
P =o J -°° 

(5.20b) + &,i,k(fc) + fcAJfefa). 

Subtracting (4.1) from (5.20) and differentiating with respect to y Ui , we get a system of integral 
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equations for e p \y) given by 





dy 



Ul 



p=0 







dy 



Ul J — oo 



F x ((y Ul ,B 2 ),t k \(y,B 2 ),t p )e^(y)dy 



(5.21a) 



+ 



8 r B * 



dllui J —OO 





F x ((y Ul ,B 2 ),t k \(B 1 ,y),t p )e p 2 \y)dy 



dy v 



p=0 







dy v 



By 



F x ((B u y U2 ),t k \(y,B 2 ),t p )e p 1 \y)dy 



(5.21b) 



+- 



a f B 2 



F x ((B 1 ,y U2 ),t k \(B 1 ,y),t p )e p 2 \y)dy 



^Uu2 J —OO 

Rewriting (5.21) with respect to k — 1, subtracting it from (5.21) and using (4.2), we obtain 



K lAp ((y Ul ,B 2 ),(y,B 2 ))e p 1 \y)dy 



(5.22a) 



p=0 

+ / K 1Ap ((y Ul ,B 2 ),(B 1 ,y))e p 2 \y)dy 

J — oo 

d 



dy Ul 



h 



4 2) (^ 2 )"E / 1 K 2Ap ((B 1 ,y U2 ),(y,B 2 ))e p 1 \y)dy 

n J —OO 



p=0 

-B 2 



+ / K 2Ap ((B 1 ,y U2 ),(B 1 ,y))e p 2 Hy)dy 
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Using (5.3), (5.19) and condition 5.1(v), and since t k -\ = t k — ft, we have 

d 



< 



dy Ul 

ht k [ B > 8_ 

, dt 



\K 1> k,t((yu 1 ,B 2 ), (y,B 2 ))\\f( X «,T 2 ) (y,t) \dy 



t=T 



h 2 d f Bl d 

+ y 



dy Ul J-oo dt 



/j 2 



Fx{(y Ul , B 2 ), t k -x\ (y, B 2 ),t)f {x ?,T 2 ) (y, *) dy 
h 2 



t=T 



< — [ifcQi.i + Si,i] := -ai,i. 

The last inequality holds applying conditions 5.1(h) and 5.1(iv) on the first term, and condition 
5.1(v) on the second term, for a suitable positive constant Si 7 i. Similarly 



d 
dy 



\Si,i,k{h) - Si )itk -i(h)\ < — [t k Qi,i + Stj] :-- 



-oil 



■ui 



for i,l = 1,2 and suitable positive constants Sn obtained from condition 5.1(v). Then (5.22) 
becomes 



(5.23a) 



l4*Wl < 



(aii + ai2)ht k " 1 rlU 



k-i 

p=0 



K hk>p ((y Ul ,B 2 ),(y,B 2 ))eU(y) dy 



+ 



•s 2 



K 1Ap ((y Ul ,B 2 ),(B 1 ,y))ef\y) 



dy 



(5.23b) \e k 2) (y U2 )\ < 



(a 2 ,i + a 2y2 )ht k 



k-l 
p=0 



if2,fc, P ((5i, 2 / U2 ),(y,S 2 )) e W(y) 



+ 



B 2 



K2,kA(B 1 ,y U2 ),{B 1 ,y))ef\y) 



Setting 7/ = max{a^i,ay}, for Z = 1,2, we can write the system (5.23) iteratively for k > 0, 
obtaining 



o (0i 



< 7l / i i 2 + ^ 1) \K 1Ap ((y Ul ,B 2 ),(y,B 2 ))\dy + hq^ \K 1Ap ((y Ul , B 2 ), (B v 

J — oo J— oo 



y))l*/ 



< ft (71*2 + + ft^9i 2) ) := ft^ 1} ; 

|e 2 2) | < ft(72t 2 + ftei 2 gi 1) + ft^f ) ) :=ft?f, 
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3 W TT„w, t3 



where we used condition 5.1(H) to bound e 2 . Here q are suitable constants independent on h 
and r. In general 



(5.24a) 



(5.24b) 



4 2) (y« 2 )l <^ 



fc-i 



fc-i 



(2) 



p=0 
fc-1 



p=0 



72^ + M^E^+^E^ 



(2) 



p=0 



p=0 



T (2) 



Since tf. < 0, from (5.24) it follows 



<7i© + fc(^Eg5 l) + ^E(i?), 

p=0 p=0 



1,2, 



and applying eq. (7.18) in [15], we get < ^i@exp[(^\ + £1)^] an d thus 

lej^GfeJI ^/iTieexpf^+es 1 )^] ! 
|ei 2) (y U2 )| <^72©exp[(C? + ^)4] • 

The result follows noting that |ej^(y u J| = 0(/i). 



(5.25a) 
(5.25b) 



□ 



Remark 5.1. T/ie numerical approximations f< y x a ,T i ){y^p) can be rewritten such that they 
depend only on f(x^,T i ){v,'tp) an d Ki,k,p{p,b,c,d), as shown in Remark 2 in [8j. Therefore, 
conditions 5.1(i) and 5.1 (Hi) are in fact assumptions on f and K. 

6. Joint distribution of (Ti, T2) for a bivariate Wiener process. Consider a bivariate 
Wiener process X solving (2.1) with constant drift /i(X(i)) = (p\, p 2 ) £ ^ 2 and positive-definite 
covariance matrix 

0-1 



pa 2 oi\J\ — p* 

for p G (— 1, 1). Then X is a bivariate Wiener process with null drift and covariance matrix 



po-\o-z 



That is, each component is marginally a Wiener process with drift pi, diffusion coefficient o~i > 
0, i = 1, 2 and p is the correlation of the bivariate Wiener process, e.g. p = corresponds to have 
independent components. For the Wiener process, the densities fx, fx? and gTi,i,j = 1, 2, % / j 
are known [9]. To apply Theorems 3.1 and 3.2, we need to calculate the unknown density f(x a ,Ti)- 
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Since is known, this corresponds to derive the conditional density fx a \Ti- The first step is 
to calculate the unknown transition density / X a , which solves the two-dimensional Kolmogorov 
forward equation 

dfx»(x,t) ^ d 2 / x »(x,t) | a 2 2 d 2 fx«(x,t) , g J d 2 / Xa (x,t) 

(6.1) _ <9/ X q(x,£) d/x°(x,t) 

with initial, boundary and absorbing conditions given by 

(6.2) lim/x»(x,t) = 5(xi - x i) 5 (x 2 - x 02 ) ; 

(6.3) lim / X a (x, t) = lim / x « (x, t) = 0; 

(6.4) /x^x,*)!^^ = /xa(x,t)| X3=Ba =0, 

respectively, where we set to = 0. The solution provided in [12] does not fulfill (6.2) when 
(Atl,/^) (0,0). Following their proof, we noted that the normalizing factor 



(6.5) exp 



(/X2/9Cri<7 2 - Ml "!) -^l + (^l/Of 1C2 - ^2Crf ) #2 

(l-p 2 )^ 2 



was missing. Since (6.5) is equal to 1 when (/zi,/i2) = (0,0), the results in [12] are correct for 
the drift less case. In presence of drift, it holds 

Lemma 6.1. The density /x« that the process never reaches the boundary B in (0,t) is given 

by 

2 

/ X o(x,t) = — — exp{K 1 (B 1 - x i) + K 2 (B 2 - x 02 )} 

(6 6) 

x exp(^ ^ t __5_j H(r|r0i ^, t)) 
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where r := r(x±, x 2 ) G (0, oo), 4> '■= <fr(xi, x 2 ) £ (0, a) and 



r 

f cos(</>) 

00 

a 



aj{B 2 - x 2 ) 2 + o\{B x - Xl y - 2a 1 a 2 p(B 1 - x x ){B 2 - x 2 ); 

a 2 (Bi - x x ) - o\p{B 2 - x 2 ), r sin(c/>) = - p 2 (B 2 - x 2 ); 

f\ ; 
0o 1 ; 

Q 2 \X\ - Q\[l 2 p 

o\o 2 {\-p 2 y 



li- 



ar ct an 



H(r,r ,(j),4> ,t) = > sin ) sin )I, 

^-^ a a 



aip 2 - G 2 ji\p 
€ (0,7r); 



K 3 = o\o 2 ^J\ - p 2 ; 



nir/a 



n=l 



K 2 t 



Here (f, (f>) are functions of (xi,x 2 ) and are obtained through a suitable change of variables 
in [12]. In (6.6) we use them instead of (x%, x 2 ) to simplify the notation. To compare our results 
with those in [12], one should introduce the transformation r = f jo\ and r$ = fo/cxi, since they 
use different constant terms. 

Remark 6.1. The distribution of the first exit time of the process from the strip (—00, B\) x 
(— co,B 2 ) is given by 

P(min(Ti,r 2 ) < t) = 1 - / 1 f 2 f^a(x,t)d Xl dx 2 , 

and it can be computed multiplying eq. (32) in [12] with the missing factor (6.5). 

Corollary 6.1. The conditional density fx a \x a ( x i>t\ x jit)> S or h3 =1,2; i 7^ j is given by 



fxf\X?{ x ii t\xj,t) 



aK^t 



-( Xj - X j)p - ( Xi - XQi) 



x exp ( —KitNj 



+ f 2 -(x j -x 0j ) 2 a 2 ( y l-p 2 ) 



2K 2 t 



(6.7) 
with 



1 — exp 



2(Bj - x j)(xj - Bj) 



a 2 t 



H(r,r ,<t>,(t>o,t), 



iVi 



ai[i 2 - a 2 p\p 
2ai 



No 



a 2 pi - o\p. 2 p 
2a 2 
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Proof. The conditional density is given by 



(6.8) 



fx?\xA%ht, \xj,t;y,s) 



/x»(x,%,s) 



and the result follows by plugging fx* given in [20] and (6.6) into (6.8). 
Then, we introduce the following 

Lemma 6.2. The conditional density fx?\Tj (%i\t) for i,j = 1,2 i ^ j is given by 



□ 



<7j7TV 2irt 



fxf \Tj{ x i\t) — 2(T3 \7cT 



X j 



exp —K, 



(Bj - x 0j )p - (xi - XQi) 



x exp —KitNj 



[po-i(Bj - x j) - crj(Bi - x i)} + °~j(B 



(6.9) 

where 

oo 

Gij{f , <j> ,xj,t) = 5., 
with 5i = 1 and 5 2 = (-l) n+1 . 



2K\t 



nirc 



Gij(r ,(j>o,Xi,t), 



rasm 



n=l 



a 



aj(Bi - xi)r 



Kit 



Proof. When Xj — > Bj, both fx a and /x a go to zero, due to the boundary condition (6.4). 
Therefore fx a \x a is indefinite, as noticed from (6.8). From the definition of <j), we have that 
(p — > ct when x\ — > B\ and <p — > when x 2 —> B 2 . In both cases, s'm(nir(J)/a) — > and thus 
H(r, fo, 0, 4>q, t) — > 0. Moreover, 



1 — exp 



2{Bj - x Qj ){ Xj - Bj) 



o, 



when Xj — > Bj. Hence, the last two terms in (6.7) produce an indefinite form. Applying l'Hopital's 
rule, we obtain 



sin 



lim 



OiOjK^JX - pH 



*i->Bi _ (■>:!}, ,:„:;,■; H:)\ 2o>(Bi - Xi )(Bj - X j) 

i exp i ^ t 



77,(5,;, 



with 5± = 1,82 = (— 1)" , depending on whether <f> — > a or (j) — >■ 0. The result follows by 
plugging this ratio into (6.7). □ 

In presence of absorbing boundaries, the joint distribution of the FPTs can be explicitly 
calculated as follows 
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Theorem 6.1. The joint distribution of (Ti, T 2 ) in presence of absorbing boundaries is given 

by 



.10) x / / / 1 1 z exp (-Kip—Xj - KiNjSi 



J-ooJsi Si^(Sj - Si) 3 

x exp \- r ° + a l} B i 2 — ~ TTYi — ~ T? J G i*( f o» ^0, Xj, s^dsjdxjdsi. 

\ ZSjxlg Z(J j\Sj Si) j I 

Proof. It follows by plugging g^ given in [20] and (6.9) into (3.1), and then simplifying the 
resulting expression. □ 

6.1. Driftless Wiener process. For a driftless Wiener process in presence of absorbing bound- 
aries, it holds 

Theorem 6.2. The joint density of(Tx,T 2 ) in presence of absorbing boundaries is given by 

.11) 



00 

X 

n=l 



^ . f nir(/) \ ( rl(tj-ti) \ 
Proof. Since pi = p 2 = 0, it follows that = K 2 = 0. Deriving (6.10) with respect to t\ 
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and t 2 , we have 



P(Ti G dh,T 2 € dt 2 ) 



E 



2tt 



2a 2 (jjtisJ{tj — t 



A3 



exp 



2t 4 Kl 



xGji(r , 0o, dxjdtidtj 



nsm 



nirc 



a 



cxp 



doc j dt^dt j 



E 



/ 27T exp 



lo w 



n sin 



a 



(6.12) 



exp 



/ fr%-^ 2 ) \ 

I I J n7r 



dhdtidtj, 



where the last equality is obtained through a change of coordinate, namely h = o~i(Bj — xj). 
The integral in (6.12) can be solved using the identity [16] 



(6.13) 



7T 



-^ h2 I u ( 1 h)dh=—e,: V 



(l 2 



setting /3 2 = (tj — Up 2 )/ (2K 2 ti(tj — ij)) ,7 = ro/(K 2 ti) and ^ = mr/a. The result follows after 
some computations. □ 

When t\ = t 2 , we have the following 

Corollary 6.2. The joint FPT density in presence of absorbing boundaries when t\ =t 2 = 
t is 



P(Ti € dt, T 2 £dt) = < 



Odt 2 
oodt 2 

(B 1 -x 01 )(B 2 -x 2) ' 



ifpe (-1,0) 
if P e (0,1) 



y%(B 1 -x 01 ) 2 +vj(B 2 -x 2)' 2 



dt 2 ifp = o 



Proof. If ti < tj, set z = tj—U, i,j = 1, 2, i / j. When 2; — >■ 0, the limit of (6.11) is indefinite, 
being of the form I„(z)/z, with v = nvr/(2a). Using the fact that I u (z) ~ (z/2) u /T(u + 1) [2], 
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we get 

if i> > I 

3.14) lim = — — -z u - 1 



h{z) = 1_ 

z-^b z 2T(i/ + 1) ' 



oo if v < 1 
\ if v = 



Since a G (0, 7r), then v > 1 for n > 2, and thus all addends in the series in (6.11) vanish for 
n > 2. For n = 2, #i = 1, 62 = —1 and thus the term in (6.11) is null, being the two densities 
symmetric. Finally, when n = 1, from (6.14), definitions of v and p, it follows that v < 1 44> a G 
(0,7r/2) ^>/)G (-1,0); ^ = 1 <^ q = 7r/2 44> p = and ^ > 1 <S> a £ (7r/2,7r) 44> p G (0,1), 
where 44> denotes if and only if. The result follows plugging the resulting expression for (6.14) 
into (6.11). □ 

Remark 6.2. To compare (6.11) with the corresponding expression in [13] for s = t\ < ti = 
t, we set 

ro 

since different transformations are used. Since 

yl — p 2 = sin a, p 2 = cos 2 a, 2(t — sp 2 ) = (t — s) + (t — s cos 2a), 

we obtain 

P(Ti G ds,T 2 Edt,s < t) 

it sin a / fn (t — s cos 2a) 

exp 1 



2a 2 v / s(t- scos 2 a)(t- s) \ 2s[(t - s) + (t - s cos 2a)] 



00 



x y(-l)"+ 1 nsin(^^)/m ( — ; uv ; " y ) dsdt. 

^ ' v a ' ^ \2s[(t-s) + (t-s cos 2a)] J 

The result differs from that in [13], that uses an incorrect identity for (6.13), as already discussed 
in [12]. 

Since the joint density f(x-,T t ) 1S n °t explicitly known, Theorem 3.2 cannot be used to de- 
termine an analytical expression of the joint distribution of (T±,T2) in presence of crossing 
boundaries. However, /(x-,t») can be numerically evaluated for bivariate Gaussian diffusion as 
described in [5]. 

7. Examples. Denote fj> the theoretical joint density of (Ti,T%) and fa its numerical 
approximation obtained applying the proposed algorithm. Here we report a brief illustration of 
the joint FPT density for two-dimensional Wiener and OU processes in presence of absorbing 
boundaries. For a bivariate Wiener, the theoretical joint density /t is given by (6.11) for the 
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Fig 1. Theoretical joint densities and contour plots of (TijT?) for two-dimensional Wiener processes in presence 
of absorbing boundaries. Parameter values common to all figures: o\ = 01 = \, p = 0.5. Chosen drifts for the top 
figures: /ii = /i2 = 0, B\ = B2 — 1 and time discretization step h = 0.005. Chosen drifts for the bottom figures: 
/ii =1,/j,2 = 1.5, -Bi = B2 — 10 and time discretization step h = 0.05. Panel (a): joint density of (T\,T2). Panel 
(b): contour plots of (Ti,Ti). 
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driftless case, and is derived from (6.10) when the drift is not null. To compare fx and fx, 
throughout we consider the mean square error (MSE), which is defined by 

MSE(/ T ) = — EE (MUM ~ MU,tj)) ■ 

i=l j=l 

7.1. Bivariate Wiener process. First, consider a symmetric bivariate Wiener process with 
null drifts, parameters a% = o~% = 1, p = 0.5 and boundaries B\ = B<i = 1. The theoretical joint 
density and its contour plot are reported in the top panels of Fig.l. The numerical approxima- 
tions are not shown, since they are indistinguishable from the theoretical. Indeed, choosing a 
space discretization step r = 0.05 and time discretization step h = 0.01 (resp. h = 0.05), we 
obtain MSE(/ T ) = 3.5859 -10" 5 (resp. MSE(/ T ) = 4.8607- 10~ 4 ). This confirms the reliability of 
the algorithm, as expected from the convergence results of the error, proved in Theorem 5.1. Not 
surprisingly, also the joint FPT density and the contour plots are symmetric, and the probability 
mass is concentrated in the area close to the diagonal t\ = t2, representing simultaneous FPTs, 
i.e. Ti = T 2 . 

Second, consider a non-symmetric bivariate Wiener process with parameters p\ = 1, pi = 2, u\ = 
oi = 1, p = 0.5 and boundaries B\ = B<± = 10. The joint density fx and the contour plot are 
reported in the bottom panels of Fig.l. They are indistinguishable from those obtained applying 
the numerical algorithm (figures not shown). As expected, the joint FPT density is not sym- 
metric and the probability mass is concentrated around the means of the FPTs and it is spread 
out according to the variance of the FPTs. Indeed, for this parameter choice, we have 

7.2. Bivariate Ornstein-Uhlenbeck process. A bivariate OU process satisfies the stochastic 
differential equation (2.1) with 

(Mi 4a j. *™=s=(:;: z 

for pi 6 M, Oij > 0, 1 < i,j < 2, ayi £ K and Yl a positive-definite matrix. 

Throughout we fix 9 = 10, o\% = l,o~i = 2, Bi = 10, for i = 1,2. First, we consider a symmetric 
OU with p\ = p2 = 1.5. The approximated joint density fx and the contour plot of (T\,T2) 
are given in the top panels of Fig. 2. With this parameter choice, the asymptotic mean pi9 of 
each component of the OU is above the boundary Bi. Also in this case, the probability mass 
is concentrated along the diagonal t\ = ti. Hence, the times when the components cross their 
boundary are similar. 

Second, we consider a non-symmetric OU with drifts p\ = 0.95 and pi = 1-5. The approximated 
joint density fx and the contour plot are reported in the bottom panels of Fig. 2. Note that the 
first component has asymptotic mean p\9 below B\, and thus the noise determines the crossings 
of the boundary. As a consequence, the probabilistic mass is concentrated in the region t\ > ti- 
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Fig 2. Approximated joint densities and contour plots of (Ti,Ta) for bivariate OU processes in presence of ab- 
sorbing boundaries. Parameter values common to all figures: <th = 022 = 2,<7i2 = 1, B% = B2 = 10. Chosen drifts 
for the top figures: i±\ = /12 = 1-5. Chosen drifts for the bottom figures: fi± — 0.95 and [ii = 1.5. Panel (a): joint 
density of (Ti,T2). Panel (b): contour plots of (T\,T2) . 

8. Conclusion. We solve the FPT problem for a two-dimensional Wiener process with 
constant drifts and non-diagonal covariance matrix in presence of absorbing boundaries. In par- 
ticular, we explicitly calculate the joint density of the FPTs and other relevant quantities, e.g. 
the transition density of the process under the boundary. For bivariate diffusion processes in 
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presence of either absorbing or crossing boundaries, explicit expressions of those densities are 
not available. Therefore, we suggest to use the proposed numerical method for the evaluation of 
f(x?,T 2 ) an d /(xj,Ti)- The algorithm does not depend on whether the boundaries are crossing or 
absorbing, and its error is shown to converge. These results are also confirmed by our numerical 
examples for the Wiener case. Using the algorithms proposed here for f(x a ,Tj) an d h 1 [5] for 
f(Xj,Ti)i h is possible to numerically evaluate the joint FPT density, provided that the bivariate 
transition density is known, e.g. for bivariate Gaussian diffusion processes. Furthermore, the 
proposed approach can be extended to the FPT problem for multivariate renewal processes with 
crossing boundaries, as argued in Remark 3.2. 

A generalization to the joint FPT distribution of a fc-dimensional process would also be of inter- 
est. However, that study requests the knowledge of the solution of a /c-dimensional Kolmogorov 
forward equation, when the process is a multivariate Wiener in presence of absorbing bound- 
aries, or of a system of k Volterra-Fredholm first kind integral equations, when the process is a 
diffusion. 

Finally, note that both the theoretical results in Section 3 and the numerical algorithm can be 
extended to the case of time dependent boundaries. 
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